load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"


;;.....................................................................................
;; nice string
;;.....................................................................................
  undef("strnice")
  function strnice(str,nc)
  begin

   string_len = sizeof(stringtochar(str)) - 1

   nstr = str

   if (string_len.lt.nc) then

   do it = 0, nc - string_len - 2
      nstr = nstr + " "
   end do

   end if

  return(nstr)
  end



begin

  ca = "AIE_1.9x2.5_EXP005"
  cb = "AIE_1.9x2.5_EXP008B"

  ya = "2006-2009"
  yb = "2006-2009"

  ya = "2009"
  yb = "2009"

  ya = "2006-2010"
  yb = "2006-2010"

  ya = "2006-2010"
  yb = "2006"

  ya = "2006"
  yb = "2006"

  seas = "ANN" 

  basepa = "~/tshr/run/run_"
  basepb = "~/tshr/run/run_"

  ma = ".cam2.h0."+ya+"."+seas+".nc"
  mb = ".cam2.h0."+yb+"."+seas+".nc"

  vname  = (/ "CLDLOW", "CLDMED", "CLDHGH", "PRECC", \
              "PRECL", "PRECT", "SWCF", "LWCF", \
              "TCF", "RFP", \
              "FSNT", "FLNT", "FSNTC", "FLNTC", \
              "FSNS", "FLNS", "FSNSC", "FLNSC", \
              "TGCLDLWP", "TGCLDIWP", "AODVIS", "AODABS", "LHFLX", "SHFLX" /)

  vprep = 24.*3600.*1000.
  vlwp  = 1.e3

  vunit = ""

  vscale = (/ 1., 1., 1., vprep, \
              vprep, vprep, 1., 1., \
              1., 1., \
              1., 1., 1., 1., \
              1., 1., 1., 1., \
              vlwp, vlwp, 1., 1., 1., 1.  /)

  vunits = "" 

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)

  print(" ")
  print(" ")
  print("..................................................................")
  print(" ")
  print(" ")
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1))
  print(" ")
  print(" ")

  fngrid = "gridarea.nc"
  fgrid = addfile(fngrid,"r")
  gridarea = fgrid->cell_area
 
  fnmask = "landmask.nc"
  fmask = addfile(fnmask,"r")
  landmask = fmask->landmask
 
;;......................................................................................
;; Open netCDF file and read in dimensions
;;......................................................................................

   filename=case(0)
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))
   
;;......................................................................................
;; read in variables & create plot
;;......................................................................................

  plot = new(Ncase+1,graphic)
 
  nv = dimsizes(vname)

  printdata = new((/4*nv+50/), "string") 
  printdata = "" 

  ip = 0 
  printdata(0) = " " 
  printdata(1) = ".................................................................. " 
  printdata(2) = " " 
  printdata(3) = " ctrl      : "+case(0) 
  printdata(4) = " test      : "+case(1) 
  ip = ip + 6 
 
  do iv=0,dimsizes(vname)-1 
  do i=0,Ncase-1
     filename=case(i)
     fi=addfile(filename,"r")

        if(vname(iv).eq."TCF") then 
             Var3dtema1 = fi->LWCF(time|:,lat|:, lon|:)
             Var3dtema2 = fi->SWCF(time|:,lat|:, lon|:)
             Var3dtema = Var3dtema1 
             Var3dtema = Var3dtema1 + Var3dtema2 
             Var3dtema@long_name = "TCF" 
          end if
  
          if(vname(iv).eq."RFP") then 
             Var3dtema1 = fi->FSNT(time|:,lat|:, lon|:)
             Var3dtema2 = fi->FLNT(time|:,lat|:, lon|:)
             Var3dtema = Var3dtema1 
             Var3dtema = Var3dtema1 - Var3dtema2
             Var3dtema@long_name = "RFP"  
          end if
  
          if(vname(iv).ne."RFP" .and. \
             vname(iv).ne."TCF" ) then 
             Var3dtema = fi->$vname(iv)$(time|:,lat|:, lon|:)
          end if   
  
          Var3dtema=Var3dtema*vscale(iv) 
  
          Var3dtemp = Var3dtema(0,:,:)
          Var3dtempland = where(landmask.eq.1,Var3dtema(0,:,:),-999) 
          Var3dtempocean = where(landmask.eq.0,Var3dtema(0,:,:),-999) 
          gridarea_land = where(landmask.eq.1,gridarea(:,:),-999) 
          gridarea_ocean = where(landmask.eq.0,gridarea(:,:),-999) 
          Var3dtemp@_FillValue = -999  
          Var3dtempland@_FillValue = -999  
          Var3dtempocean@_FillValue = -999 
          gridarea_land@_FillValue = -999  
          gridarea_ocean@_FillValue = -999 
  
          global_total = sum(Var3dtemp*gridarea)/sum(gridarea) 
          global_total_land = sum(Var3dtempland*gridarea_land)/sum(gridarea_land) 
          global_total_ocean = sum(Var3dtempocean*gridarea_ocean)/sum(gridarea_ocean)      

          if(i.eq.0) then 
             printdata(ip) = strnice(vname(iv),15)+" : "+sprintf("%9.4f", global_total) +"   "+sprintf("%9.4f", global_total_land)+"   "+sprintf("%9.4f", global_total_ocean) 
             print(""+printdata(ip))  
             ip = ip + 1 
          end if 

          if(i.eq.1) then 
             printdata(ip) = strnice(vname(iv),15)+" : "+sprintf("%9.4f", global_total) +"   "+sprintf("%9.4f", global_total_land)+"   "+sprintf("%9.4f", global_total_ocean) 
             print(""+printdata(ip))  
             ip = ip + 1 
          end if 

     if(i.eq.0) then
        vara = Var3dtemp
        varab = Var3dtempland
        varac = Var3dtempocean 
     end if

     if(i.eq.1) then
        varb = Var3dtemp
        varbb = Var3dtempland 
        varbc = Var3dtempocean 
     end if

     delete(Var3dtemp)

  end do

  varc = vara
  varc = varb - vara

  vard = vara
  vard = varbb - varab

  vare = vara
  vare = varbc - varac

  global_total = sum(varc*gridarea)/sum(gridarea)
  global_total_land = sum(vard*gridarea_land)/sum(gridarea_land)
  global_total_ocean = sum(vare*gridarea_ocean)/sum(gridarea_ocean)

  printdata(ip) = strnice(vname(iv)+" diff ",15)+" : "+sprintf("%9.4f", global_total) +"   "+sprintf("%9.4f", global_total_land)+"   "+sprintf("%9.4f", global_total_ocean) 
  print(""+printdata(ip)) 
  ip = ip + 1 

  printdata(ip) = "  " 
  print(""+printdata(ip)) 
  ip = ip + 1 


  end do
  
  asciiwrite("table_"+cb+"_"+ca+"_"+ya+"_"+yb+"_"+seas, printdata) 

end


